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foe EN ERODUCTION 


The need for readily computed parameter estimates is 
great. Maximum likelihood estimators are known to be asymp- 
totically efficient, but in many Settings they are hard to 
find. The most popular alternative is the method of moments 
which usually yields readily computed estimates, but some- 
times these estimates are not very efficient. See ref. [2]. 
This study looks at the efficiency of a method which uses 
the probability generating function (empirical La Place trans- 
form evaluated at a desirable value of its argument, uo): 

The method presented herein requires computing power 
greater than that of the method of moments, but less than 
that of maximum likelihood, which calls for a psi function 
capability. It can be used with a hand held calculator that 
has the lower order transcendental functions, l1.e., logarithm, 
roots and power. 

The basic idea 1S to select a system of eStimating equa- 
tions which equates various Eeablet les to their expected 
values. In chapter II, estimating equations were set and 
efficiency computations were done using theoretical work’ 
presented in [3]. Chapter III contains the procedure to be 
followed, applications of the method to some data sets, and 
comparisons with the other methods. 

In chapter IV, efficiency and optimum uy and to (=-log u,) 
tables for various p-r combinations are given and efficiency 


contours in three different planes are graphed. 





II. METHODOLOGY 


A negative binomial random variable X has the probability 


function 


Pe x) ox 


£ (X=x;3r,p) ees 


Zee) 


POR w—"O0,l,2,0.06 0 < b, 0 < p< 1, p+ q= 1. The partial 


derivatives of its logarithm are, using h = log f 


oh - ©£_x 
dp Pp q 
ae == Cet ee 2 (ro slog p 


Using the basic recursive formula for the psi function [1] 
ee aay 


< 
Wir+x) - ¥(r) = } 


one may then express the system of maximum likelihood 


equations as 


x ~ ae | = @) 
Pp 
(24.3) 
Ave ) t teleg pp == 0 
ee i a Were | 
od eae e fl 





where the values x x are the data that result when 


loots 
the population (2.1) is sampled n times. The system (2.3) 
Monon linear In f,o and Gilf£rrcult to solve. iterative 
methods must be used and the second member of (2.3) (or 
some version thereof) must be recomputed in each cycle. 
This is the main computational difficulty in uSing maximum 


likelihood in this Setting. 


The information matrix A, which is defined as: 





3°h 3“h 

ee dpar 
A = -E 

3“ 3“ 

drdp yr? 


J 


1s needed for efficiency calculations and can be developed 


using the methods presented in {2]. The result is 
pcs es 
Gp a 
A = 
IE 
D Noo 
where 
pees 1 az ; 
An = W'(r) Be 2 be) ) 





The determinant of A is found (using [2]) as 





1 es Gia: oe ga 
VA 5 3 hen ° ers Seley 
i n=1 


The two estimating equations or and Jor which are 


exploited herein, are given as 


g, (x-prX) = Sees = 0 
(2.5) 
ees. We Pant 
J, (L,P7x) le (T=9t? 0 
where 
Come 2 es 
n 
1=1 
is the empirical generating function. (Note: the substitution 
t = exp(-u) converts this to La Place transform). Reduction 


to a Single equation is effected by solving 9, = 0 for p 


in terms of x and ee 





p = (2.6) 


and substituting this into J5- A function f(r) is obtained 


x — 
in terms of known variables, namely t , x, and t, 
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2G) = ee a een (2.7) 


r+ x(l-t) 


* 


We can solve equation (2.7) for r by Newton's method, the 


iterative relationship being 


where 


re ae - 
r = p(—xU-t) + log —— ee 


Sat) ere ext Es 
jg awe (1s Pacts (lb =o) r+ x(1l-t) 


assuming an initial value of r and a suitable t can be 
found (see table 4.2). After solving for r, p can be found 
Wsang (226): 


Efficiency for this r-p pair is given by (see [2,3]) 


a) 
ere = IM 1: (2.8): 


A 


where 


2 
im} = Jal” (229) 


jC | 


|A| is given in (2.4), the matices A and C in (2.9) are 


defined as 


AAR 





OG, 9D 4 





dp or 
A = &£E: 
99 5 95 5 
op or 
Fy 3199 
C = neE 
2 
3139 39 
and J11G5 are from (2.5). After some calculation, A and 
CG Ellrn Out EO be 
=. = 
2 
D Pp 
A = (2.10) 
= = = 
-(7P2* 1 ag os (GE e) * tog (FS) 
(l=qe) 
rg _ rqp™"(1-t) 
+ 
p (1-qt) **1 
Gone (2 410) 
r-1 
= eel ee (a ee ee 
Mee” Teaes fF he 


The efficiency of the estimation scheme is found, substituting 


[A], |C| and JA] into (2.9) and (2.8), as 


i, 





ee Here eel ven ee 
ie) ISG 





2) Ceti) 


zi co n ' 
| (lat) ox on ee pee Cy eet 
[(——+—5— -p (1-qt) “- rqp (1-t) th (m#1) (xt+n)? 


(ZZ) 


Detailed computations are given in Appendix A. The initial 
value of r is found uSing the method of moments (D,r), 
and the value of t used maximizes the efficiency (2.12) 


when evaluated at D, r. 
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IIIf. USAGE AND APPLICATION 


In this chapter the procedure to be followed, in order 
to use the method developed, iS given in greater detail and 
is applied to three sets of data taken from [4] and sixty 
shee of data generated by computer simulation uSing variqus 


combinations of the parameters. 


A. PROCEDURE 
The steps of procedure are given below. 


1. Using the method of moments, find D and r as 





starting values. Starting values are given by (3.1) 
eseemih2)))) 
Dp = x/s7 
(331) 
* <2 
py = 
s*-x 


where x is the sample mean and 5? is the sample variance. 

2. Using D, r found in step 1, findthe value of t, 
which maximizes efficiency function (2.12), and call it 
toe Brabealis6 find the efficiency, which iS maximum with this 
pair D, r and toe and call it EFF1. Tables (4.1) and (4.2) 
have been prepared for this. Single variable search and 
golden section search were used by the author to find ee 


and EFF1. 
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Bee USIANG rasa Starting value and toe found in step 2, 
find r* as a solution of (2.7), and p* uSing (2.6). Newton's 
method was used by the author to find r*. 

4. Step 2 can be repeated with p* and r* replacing p 
and r, and using a new ty (update to from Table (4.2) if 
the current estimated efficiency (Table 4.1) is not suffi- 


ciently high). This step involves the recomputations of 


t*~ and is seldom needed. The new SEETeELency USing p~7 727" 


and updated te is called EFF2. 


5. Stop if all the following conditions are met. 


Ga Phe) eseenes fer |EER2 EFFI) = ¢> 


and 


~s 


Ip* - p| 


{A 
'¢)) 


ee 
and 


ro 5 


[A 


fa 


Then p* and r* are the estimators of true parameters p and 
r, and the estimated efficiency is EFF2. 
6. If one or more of the conditions in Step 5 are not 


met, let 
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EFFL = EFF2 


and go to step 3 and repeat the steps following until the 
conditions in step 5 are met. Cy, @57 2, and ey are the 


stopping criteria which can be chosen by the uSer. 


B. EXAMPLES 

Examples one through three are the applications of the 
method to the Cricket score data of Reep, Polard and 
Benjamin [4]. Example four 1s based on data generated by 
computer and contains sixty cases. 

Example 1: Applying method developed to the Cowdrey 


data [4] the following information 1s obtained. 


Sample mean = 1.692 
Sample variance = 4.343 


O32 


Joe 


r* 0.816 


Optimum ty 0.459 


0.996 


Efficiency 


Number of Iterations = 3 


With the same data maximum likelihood estimators D, r, 


and method of moment estimates D, r are given in [2] as 


16 





Os 31 


FY > 
ll 


Oh S78 


O 
ll 


0.390 1.080 


O 
lI 
ry 
lI 


Example 2: Example 1 was repeated for Barrington data 


{4} and following information 1s obtained. 
Sample mean = 2.095 
Sample variance = 4.939 
p* = 0.346 
Cae Leet 


- 938 


ll 
ro) 


Optimum te 


Efficiency 0.995 


ll 
WwW 


Number of Iterations 


0.345 1.014 


Ky 
ll 


P 


0.424 ig 1.543 


TH 2 
ll 


Example 3: Like the previous two examples, the following 
information is obtained by applying the method to the 


Graveney data [4]: 


Sample mean = 1.570 

Sample variance = 4.474 
p2 = 0.315 
r*¥ = Q.722 


Ry 





Optimum te 0.430 


Efficiency = 0.996 

Number of iterations = 2 
= = 0.317 r = 0.729 
p = 0.351 r = 0.849 


When the first three examples are examined carefully 
it 1S noticed that the method developed in this study is 
almost as efficient as maximum likelihood method. Estimators 
are pretty close to maximum likelihood estimators and 
much better than those found by the method of moments. 
Example 4: This example is based on the data generated 
by computer. Sixty cases, each of which has a different 
sample size and parameters p and r. Four sample sizes, 
which are 15, 30, 50 and 100; three r values; 0.5, 2.5 and 
5; and five p values, 0.05, 0.1, 0.3, 0.5, 0.8 are used. 
In each case the data 1s generated for given sample size 
(n), p and r. Results, which contain case number (case no.), 
Sample size (n), p, r, method of moments estimates (p,x), 
ped ram Uy at p-r (Initial us), efficiency at the end of 
the first iteration (Initial eff.), empirical La Place 


transform estimates (p*,r*), optimum u_ at p¥*-r* (Final Uo). 


O 
efficiency at p*-r*-final Uy (final eff) and number of 
iterations (No. of itr.), are tabulated in the following 


pages. Data generated for the sixty caSeS are given in 
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Appendix B. The computer program, written in FORTRAN IV, 
is also given in Appendix C. 

When all the cases are examined, it's Seen that esti- 
mates are not given for cases 2, 3 and 18. Case 3 provides 
nothing to work with as all samples are zero, yielding 
X = 0, s° = 0. The other two cases have s¢ < x which 
Signals trouble because the method of moments estimates 
: < 0, p > 1; values which cannot be used to initialize 
our iteration scheme (3.1) or the maximum likelihood scheme 
(223). 

For each of these cases the function f(r) of (2.7) 
is a non-negative decreasing function of r which tends to 
+o aS r + 0 and is asymptotically zero as r > ~. Thus no 


root exists. The maximum likelihood approach yields a 


comparable situation. 


ee, 
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IV. RESULT AND CONCLUSION 


cre tt bass chapter efficiencies for various p-r pairs were 
given in Table (4.1). The method developed was used to 
Peiicuiate table entries. Three entries were given for each 
p-r pair. Entry in the middle is efficiency value for that 
pair at optimum t, or ul: Entries inside parenthesis are 


also efficiencies and found as follows: 


Let na log te 


If we miss Us by +50 percent, what would the efficiency be? 


Efficiency inside the parenthesis above the optimum efficiency 


answers the question for u, = 0.5 Use and other one below 
the optimum answers for uy = aes Uo- 
The change u = -log t (or rather its inverse t = e ”) 


changes the probability generating t” function into the 
empirical La Place transform. It is seen that a 50 percent 
error in its argument u results in only a minor degradation 
of efficiency, whereas a 50 percent error in t would have 
disasterous effects. In all of the examples, the initial 
u, was well within this 50 percent range, the worse case 
being 38 percent. 

After examining Table (4.1) one can conclude that effi- 
ciency is monotonically increasing with increasing p and 


increasing r for r greater than one and minimum efficiency 


at optimum t, or uy 1s greater than 0.95. 
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Following Table (4.1) optimum to and Uy (in parenthesis) 
were tabulated in Table (4.2). One may observe that t, 
decreases with p and increases with r, whereas un behaves 
oppositely. 

After Table (4.2) efficiency contours were graphed in 
p-r plane (Figure 4.1), in mean-r plane (Figure 4.2) and 
in mean-exp(r) plane (Figure 4.3). Compared with the 
Figure in [5] Figure (4.3) always gives better efficiency 


values than all the other alternatives. 
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Figure 4.3 
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Efficiency contours in MEAN-EXP(r) plane 
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APPENDIX A 


DETAILED COMPUTATIONS 


The two estimating equations Jy and J given in chapter 


II are 


-x-id = 
oi ee 0 


x Daee 


Go =t — GEge’ = 0 


Matrix A's entries were given as follows: 








ee. ego : 1 = 2 
a= B(d) = (ag) = dl 
OS Be . Se mer 
: le esas} 
(l=-qt) 
A,5 = Be) = E(in (22a =—a=) (ye el | 
= in(—at =a) (7 E)* 
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Matrix C's 


11 


12 


entries were given as follows: 





nE(S}) = n EL(x - 23%) 
n Var(x) = on Var (z= ae.) 
n i 
i=] 
n 
ny Var( } en) eee eile (cs) 
n 
1=1 
rq 
2 
Pp 
Be c\\ ee 
n E[ (x ay) Os (7-g¢? ) J 
ee n ar 
n Cov(X,t ) = on Cov(= i’t ) 
i=1i 
Cov ( x. ,t*) 
n n n 
x 1 oo 
} Cov(x,,t™) = ) Cov(x,,=> J t -) 
i=l i=l j=1 
n on 
1 x. 
= ) ) Cov(x,,t J) 
i=l j=l 
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x 
Since Cov (x,,t J) = 0 for 1 # j because of independence, 


it follows that 


>< 
mk 


n 
- = 
Cc =e = ) Cov(x,,t ~) 


12 Zt 
1=1 


= Cov(x,t*) = E(x-t*) - E(x)E(t*) 


One can calculate E(x-t*) uSing the following relation. 


G(t) = E(t*) G'(t) = E(xt*7?) 
Then 
E(xt™) = tG'(t) 
where 
Cie) 0S (Gees) 


taking the derivative and multiplying by t, E(xt*) is found 


as 


trqp” 


E(xt*) = 
(l-qt) tt 
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Then it yields 











C Eye 7 aeiow ue Lae 
12 21 Heat) 7 te 
7 2: ax _ Dear. 2 
C5 = n E(g5) = E[(t (T-ge" | 
n 
xX i et 
= nVar(t') = on Var (— y ce a) 
1=1 
n 
er 
= =var( ) t 7) 
1=1 
= Var(t*) = B(t?%) = fee) 1 
Secs = ce 1 
2x 
= (—25)" - (Sp 
eee L-qt 


As the reader will recall, |A| was given by (2.4) and can be 


written as follows. 


Call the summation from n = m+tl to m the remainder. The 
remainder can be adjusted such that it will be smaller than 


any specified small value (e). 
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In! 
Because nat < 1 (reciprocal of binomial coefficient), 
. m+l m+1 
y 2 i < g ) ope = eee < > 
n+l — m+2 (m+2)p — 
n=m+1 “n=0 


A value of m can be found such that 


EL 
qm 


(mt+2)p — © 


holds. Then this m guarantees the remainder will be smaller 


than ¢«, and the efficiency equation of (2.12) becomes 


pizqt in (274k) 7 qg(i-t))2 





q r+ m on ; 
(1-qt) 2 xr 2 q rin! 
[ -p (l-qt)“-rqp (1-t)“]( } = ) 
Givces)> nap ntl (r+n)! 
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APPENDIX B 


SAMPLES GENERATED BY COMPUTER 


Samples generated by computer for sixty cases used in 
chapter three, example four, are given in the following 
pages case by case. 

See subroutine GNRNB of computer program-3 of Appendix 


C for the simulation procedure, which was taken from [6]. 
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CASE 4 
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CASE 7 
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14 1 
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19 i 
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CASE 8 
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B 3 
5 iL 
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ad 1 
8 3 
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ala a 
14 1 
. CASE 9 
0 5 
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4 2 
5 2 
6 2 

CASE 10 
HAE i 
14 1 
74a il 
30 1 
33 2 
39 a 
40 1 
57 iL 
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VAL. PRO 
CASE 10 
(Cont'd) 
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CASE 1l 
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CASE 30 
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VAL. 
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bee 
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CASE 3l 


CASE 32 


CASE 33 


CASE 34 


‘¢ 
ne 
WY 
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WW 
U FPP ARN ARONW OHO 


FP FN FO U1 OO OO OO 
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VAL. 
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CASE 36 
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CASE 38 
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se 8)e 


VAL. 
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re 
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CASE 39 


Im Hm Wo 0 WM ~!) 


CASE 40 


ee Oe ee Oe ee Oe Oo OO ee Ce elt Oe ol all ool oe 
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BPRRPN PRP RPP POPE PEP EPNi PONE BRE RPE WWWne 


QQ 
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~/ 
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i 
NO 
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ee ee ee ee ee oe oe On ee ee ee Oe ee 


al 


“SIO & WN FO 


1 & WN FO 


m WN fF © 


NFO 


CASE 46 


CASE 47 


CASE 48 





VAL. FRQ 
CASE 49 

0 1 

1 8 

2 12 

3 13 

4 5 

fs) 1 ay 

6 11 

is 9 

8 4 

9 2 

10 2 
a1 fs) 
ils, 7s 
20 1 
CASE 50 

0 Bigs 

a 21 

2 22 

3 17 

4 12 

is) 5 

6 Z 

7 2 

8 3 

is 1 
CASE 51 

0 70 

il 18 

2 8 

3 1 

4 2 

6 1 

CASE 52 

0 a! 

2 3 

3 3 

4 4 

S) 4 

6 3 

qu 7 

8 eZ 

9 y 


VAL. FRO. VAL. FRO. 
CASE 52 CASE 54 
(Cont'd) 
0 60 
10 7 ik 15 
11 2 2 10 
12 5 3 4 
13 5 4 6 
14 4 5 2 
is é 6 1 
1G 5 3 i 
Gi 4 14 1 
18 1 
19 1 CASE 55 
20 9 
at t 13 1 
22 2 
24 2 15 3 
18 it 
ae 19 1 
27 1 
29 1 Z2 1 
23 1 
24 2 
CASE 33 35 4 
26 2 
0 5 28 1 
7 9 29 i 
2 3 30 i 
3 10 21iL 1 
4 1G 32 2 
5 9 33 2 
6 7 34 1 
7 7 36 3 
8 7 BG 6 
9 3 a7 1 
10 7 38 4 
11 5 39 3 
12 2 40 2 
ig 2 41 1 
14 2 42 2 
is 2 43 1 
a7 1 44 2 
19 1 45 iL 
23 7 46 3 
24 1 47 2 
49 1 
50 2 
51 1 
52 1 
53 3 
54 2 
BS 2 


SZ 


VAL. 


CASE 55 
(Cont'd) 


CASE 56 


PNHUPWEHRPNN WAN PNYWN BN WN H~ NE 


ERO. 


SN ee Ne ne nn el 





VAL. FRQ. VAL. BRO. VAL. BRO; VAL. ERO, 


CASE 56 CASE 58 CASE 58 CASE 59 
mca) a miGoneia) (Ganesh 
21 1 
27 4 33° 2 116 1 30 3 
28 2 35 1 118 , 32 3 
29 2 37 1 119 1 34 2 
30 2 43 1 121 1 35 5 
31 3 44 2 122 1 36 1 
33 3 45 4 124 1 37 2 
34 3 46 4 125 1 38 1 
35 1 50 1 127 2 39 1 
36 1 53 1 131 1 40 2 
37 1 54 1 135 1 Al 2 
38 1 55 2 137 1 42 2 
39 2 57 1 140 1 A4 1 
40 2 60 1 141 1 45 1 
46 1 64 1 144 1 47 1 
48 1 65 1 147 1 49 2 
50 i 68 1 152 1 51 1 
53 1 69 3 153 1 52 1 
54 1 71 i ie 1 53 1 
72 1 72 2 180 1 54 1 
16 il 73 1 185 2 56 1 
74 2 190 1 57 2 
CASE 57 16 2 192 1 58 1 
are 77 2 193 1 61 1 
78 1 197 1 63 1 
: 3 79 1 64 1 
; He 80 1 CASE 59 65 1 
: : 81 1 66 1 
82 1 67 2 
4 9 7 2 
33 1 69 1 
5 7 10 1 
85 2 70 il 
6 2 11 2 
88 2 72 1 
7 2 12 1 
90 2 73 1 
g 4 13 1 
91 1 73 2 
9 1 Ls 1 
93 1 81 1 
10 2 17 1 
95 1 82 3 
12 2 18 1 
96 2 85 2 
13 1 19 1 
97 1 8g 1 
14 1 20 2 
99 1 90 1 
17 3 21 2 
100 3 91 1 
18 3 22 4 
101 3 93 2 
21 1 23 3 
103 1 98 1 
23 1 24 2 
104 1 102 1 
27 1 25 2 
105 1 103 1 
29 1 26 1 
ag 1 108 2 27 2 oe i 
109 2 28 1 110 1k 
110 1 a : 140 1 
114 1 
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VAL. 


Oo OND UI S&S WN FF’ © 


CASE 60 


PRE RPRPNONFPWHENNNNNKE BHPNRFPRP UPB NWHABRUD 


FRQ. 
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APPENDIX C 


COMPUTER PROGRAMS 


Computer programs, to prepare Table (4.1) and Table (4.2) 
(Computer program 1), to find the estimates for the examples 
one through three of chatper three (computer program 2), 
and to generate the data for the sixty cases of the example 
four of chapter three and to find the estimates for them 
(computer program 3), are given in the following pages. 

There are two subroutines, namely GOLDEN (Golden section 
search), SVS (Single variable search), and one function, 
namely FCN (to evaluate efficiency function), which are 
common to the three programs mentioned above. TheSe are 
given at the beginning under the title COMMON SUBROUTINES. 
After that the three programs are given. 

One who wants to use any one of the programs, must include 
the common subroutines after the main program of the computer 


program he chooses. 
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gS 


aly 
18 


Pai 
999 


COMMON SUBROUTINES 


SUBROUTINE SVS(FCN,TNOT) 
IMPLICIT REAL*8 (A-H,O-Z) 
COMMON P,R,DETLAM 
Q=1.O0DO-P 
UPPER=1. 0DO/DABS (Q) 
ALOW=-1.ODO/DSQRT (DABS (Q) ) 
RIOU=0.000001D0 
START=0.9999 
U=START 
S=0.0001D0 
A=FCN (U) 
U1L=U+S 
B=FCN (UL) 
IF (B.GT.A) GO TO 5 
IF (B.EQ.A) GO TO 6 
M=3 
U2=START+S*2.QD0** (M—2) 
IF (U2.GE.UPPER) GO TO 99 
A=B 
B=FCN (U2) 
Peano.) GO TO 7 
IF (A.LT.B) GO TO 8 
M=M+1 
U=U1 
U1L=U2 

GO-TO vi 
CALL GOLDEN (FCN,UL,U2,RIOU,50,TNOT) 
GO TO 999 


CALL GOLDEN (FCN,UL,UPPER,RIOU,50,TNOT) 


GOyEOn I 99 
CALL GOLDEN (FCN,U,U2,RIOU,50,TNOT) 
SO7l© 997 
M=2 
CALL GOLDEN (FCN,U,U1,RIOU,50,TNOT) 
GO TO 999 
M=3 
U2=START-S*2.0D0** (M-3) 
IF (U2.LE.ALOW) GO TO 77 
B=A 
A=FCN (U2) 
IF (A.EQ.B) GO TO 17 
IF (A.GT.B) GO TO 18 
U1L=U 
U=U2 
=M+1 
GO TO 2 
CALL GOLDEN (FCN,U2,U,RIOU,50,TNOT) 
GO-TO 999 
CALL GOLDEN (FCN,U2,U1,RIOU,50,TNOT) 
GOeLO= 292 
CALL GOLDEN (FCN,ALOW,U,RIOU,50,TNOT) 
RETURN 
END 
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20 


30 
40 


SUBROUTINE GOLDEN (FCN, ENDL, ENDR,RIOU,NMAX, XOLD) 
IMPLICIT REAL*8 (A-H,O-Z) 
COMMON P,R,DETLAM 
SIGMA=(3.0D0-DSQRT(5.ODO) )/2-.0D0 
AL=ENDR-ENDL 
XOLD=ENDL+SIGMA*AL 
FOLD=FCN (XOLD) 

DO 1 I=2,NMAX 
XNEW=ENDL+ENDR-XOLD 
FNEW=FCN (XNEW) | 
IF (FNEW.LE.FOLD) GO TO 5 
IF (XNEW.LT.XOLD) GO TO 6 
ENDR=XNEW 

AL=ENDR-ENDL 

GO TO 30 

ENDL=XNEW 

GO TO 20 

IF (XNEW.LT.XOLD) GO TO 7 
ENDL=XOLD 

XOLD=XNEW 

FOLD=FNEW 

GO TO 20 

ENDR=XOLD 

GO TO 8 

IF (AL.LE.RIOU) GO TO 40 
CONTINUE 

RETURN 

END 


FUNCTION FCN (X) 

IMPLICIT REAL *8 (A-H,O-2Z) 

COMMON P,R,DETLAM 

ZR =1.0D-73 

Q=1.0D0-P 

A=R*P** (R+2.D0) /Q 

B=1.0D0-Q*X 

C= (DABS ( (B*DLOG (B/P) /P—Q* (1.0D0-X) /P) )) **2.D0 
D=B** (2.D0*R+2.0D0) / (1.0D0-Q*DABS (X) **2.D0) **R 
hPa bi 2. DOTk Or ~R~ (DABS (]. OD0=—xX) )**2.D0 
=-A*C/DETLAM 

F=D-E 

IF (DABS (F).LT.ZR) GO TO 5 

FCN=G/F 

GO TO 7 

FCN=-1.0D0 

RETURN 

END 
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COMPUTER PROGRAM 1 


IMPLICIT REAL*8 (A-H,O-Z) 

EXTERNAL FCN 

DIMENSION PR(14) ,RR(26) , EFF (26,14) , TZERO(26,14) , TTEST (52) ,EFF1(26, 
114) ,EFF2 (26,14) 

COMMON P,R,DETLAM 

DATA EFF’, TZERO/728*0.0D0/ 

DATA EFF1,EFF2/728*0.0D0/ 

DATA PR/.05D0, .06D0, .07D0, .08D0, .09D0, .1D0, .2D0, .3D0, .4D0, .5D0, .6D 
10,.7D0,.8D0, -9D0/ 

DATA RR/.5D0,1.D0,1.5D0,2.D0,2.5D0,3.D0,3.5D0,4.D0,4.5D0,5.D0,6.D0 
1,7.D0,8.D0,9.D0,10.D0,11.D0,12.D0,13.D0,14.D0,15.D0,16.D0,17.D0,18 
1.D0,19.D0,0.0D0,0.0D0/ 

DATA TTEST/52*0.0D0/ 

WRITE (6,444) 

444 FORMAT ('1') 

EPS=0. O0O001DO 

DO 1 J=1,14 

P=PR (J) 

O=1.0D0-P 

K=1 

7 TEST=0** (K+1) /(@ * (DFLOAT (K) +2. 0DO) ) 

K=K+1 

GO TO 7 

8 M=K 
WRITE (6,21) M 
21 FORMAT (15X, 'UPPER BOUND=' , 18) 

DO 2 I=1,24 

R=RR (I) 

DETLAM=0 . ODO 

DO 10 N=L,M 

RATION=1 . ODO 

DO 12 JJ=1,N 
12 RATION=RATION*DFLOAT ) / (DFLOAT (JJ) +R) 
DETLAM=DETLAM+Q* *N*RATION/DFLOAT (N+1) 
10 CONTINUE 

CALL SVS (FCN, TNOT) 

TZERO (I ,J)=INOT 

TPL=DEXP (0. 5D0*DLOG (TNOT) ) 

TMIN=DEXP (1. 5DO*DLOG (‘INOT) ) 


EFF 2 (I,J)=-FCN (TPL) 
EFF 1 (I,J) =-FCN (TMIN) 
EFF (I,J) =-FCN (TNOT) 
2 CONTINUE 
1 CONTINUE 


WRITE (6,200) (PR(I) ,I=1,14) 
200 FORMAT (///,50X,'P VALUES',/,12X,14 (2X,F5.3,1X) ,/,12X,14 (2X, '--—--——- 
1’ ,1X),/) 
DO 13 I=1,24 
WRITE (6,202) (EFF1(I,J) ,J=1,14) 
WRITE (6,201) RR(I), (EFF (I,J) ,J=1,14) 
WRITE (6,202) (EFF2(I,J) ,J=1,14) 
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201 


202 
203 
13 


FORMAT (3X, 'R=',F5.2,2X,14 (2X,F5.3,1X) ) 
WRITE (6,203) (TZERO(I,J) ,J=1,14) 
FORMAT (12X,14(1X,'(',F5.3,')')) 

FORMAT (5X, 'OP. TNOT', 14 (2X,F5.3,1X) ,//) 
CONTINUE 

STOP 

END 
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COMPUTER PROGRAM 2 


IMPLICIT REAL*8 (A-H,O-2) 
EXTERNAL FCN 

DIMENSION FRQ(20) ,OBS (20) 
COMMON P,R,DETLAM 


WRITE (6,444) 
444 FORMAT ('1') 
EPS=0 .0001D0 
WRITE (6,91) 
ei FORMAT (5X, 'READ IN NO. OF OBS. AND VARIATES IN 215 FORMAT',/) 


READ (5,81) NOBS,.NVAR 
81 FORMAT (215) 

WRITE (6,92) 
92 FORMAT (5X, 'NOW GIVE ME VARIATE VALUES AND CORRESPONDING FREQUENCY 

1IN THE FORMAT 2F10.5 FOR EACH CARD',/) 

SUM=0 . ODO 

SUM1=0. ODO 

DO 4 I=1,NVAR 

READ (5,82) CBS (I) ,FRQ(T) 
82 FORMAT (2F10.5) 

SUM=SUM+OBS (I) *FRQ (TI) 

SUM1=SUM1+ (OBS (I) **2.0D0) *FRQ(T) 

XBAR=SUM/DFLOAT (NOBS) 

SSQR= (SUM1-DFLOAT (NOBS) *XBAR**2 .0D0) /DFLOAT (NOBS~1) 

WRITE (6,209) XBAR, SSOR 
209 FORMAT(15X,'USING GIVEN DATA FOLLOWING INFORMATION IS OBTAINED' ,// 

1,20X,'SAMPLE MEAN =',F7.4,//,20X,'SAMPLE VARIANCE =',F7.4,// 

2) 

P=XBAR/SSOR 

R=XBAR**2.0D0/ (SSQR-XBAR) 

_Q=1.0D0-P 

K=1 
7 TEST=0** (K+1) / (DFLOAT (K+2) ) 

IF (TEST.LE.EPS) GO TO 8 

K=K+1 


DO 12 JJ=1,N 
12 RATION=RATION*DFLOAT (JJ) / (DFLOAT (JJ) +R) 
DETLAM=DETLAM+Q**N*RATION/DFLOAT (N+1) 
10 CONTINUE 
[=] 
95 CALL SVS (FCN, TNOT) 
EFF 1=-FCN (TNOT) 
RSTAR=R : 
PSTAR=P 
SUM2=0 . ODO 
DO 14 J=1,NVAR 
14 SUM2=SUM2+F RQ (J) * INOT**OBS (J) 
TXBAR=SUM2/DFLOAT (NOBS) 
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77 


88 


52 
51 


29 


999 


30 


66 


17.4,//,20X,aFFICIENCY 


CALL COMPR (TNOT, TXBAR, XBAR, RSTAR) 
P=R/ (R+XBAR) 

Q=1.0D0-P 

K=1 

TEST=0%* * (K+1) / (P*DFLOAT (K+2) ) 

IF (TEST.LE.EPS) GO TO 88 


K=K+1 

GO TO 77 

M=K 
DETLAM=0.0DO0 | 
DO 51 NE1,M 
RATION=1 . ODO 
DO 52 JJ=1,N 


RATION=RATION*DFLOAT (JJ) / (DFLOAT (JJ) +R) 
DETLAM=DETLAM+Q**N*RATION/DFLOAT (N+1) 
EFF 2=-FCN (TNOT) 

TF (DABS (P—PSTAR) .GT.0.005) GO TO 29 

TF (DABS (R-RSTAR) .GT.0.005) GO TO 29 

IF (DABS (EFF1-EFF2) .GT.0.04) GO TO 29 
WRITE (6,94) P,R,EFF2,TNOT,L 
FORMAT (20X, 'ESTIMATE OF P ' F7.5,//,20X, ESTIMATE OF R 
' F7.5,//,20X, ‘OPTIMUM TNOT 
27.5,//,20X, 'NO.OF ITERATIONS =',F14,//) 

GO TO 999 

L=L+1 

GO TO 95 

STOP 

END 


SUBROUTINE COMPR(T,TX,XB,RST) 
IMPLICIT REAL*8 (A-H,O-Z) 
COMMON P,R,DETLAM 
R1=RST 
B+R1+XB* (1. 0D0-T) 
TY=(RL/B) **RL 
F=TX-TY 
FP1=-TY 
FP2= (XB* (1. ODO-T) /B) +DLOG (RL/B) 
FP=FPL*FP2 © 
R2=R1-F/FP 
IF (DABS (R2-RL) .LE.0.00001D0) GO TO 66 
R1=R2 
GO TO 30 
R=R2 
RETURN 
END 
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541.. 


978 


979 


444 


S21 


390 


370 


380 


COMPUTER PROGRAM 


IMPLICIT REAL*8 (A-H,O-Z) 


EXTERNAL F'CN 


DIMENSION FRQ(200) ,OBS (200) ,PVECT (10) , RVECT (10) 


1,1SD(200) 
COMMON P,R,DETLAM 


DATA RVECT/5.0D0,2.5D0,0.5D0,7*0.0D0/ 
DATA PVECT/.8D0,.5D0,.3D0, .1D0, .05D0,5*0.0D0/ 


DATA NVECT/15,30,50,100,6*0/ 


CALL OVFLOW 
WRITE (5,541) 


3 


|, NVECT (10) , DAT (200) 


FORMAT (//,15X,'READ ISEED USED IN GENERATING DATA’) 


READ (5,978) IS 
FORMAT (I10) 
WRITE (6,979) I5 


FORMAT (15X,'ISEED USED TO GENERATE DATA=' ,110) 


CALL INT (IS,ISD, 200) 


WRITE (6,444) 
FORMAT ('1') 
EPS=0.0001D0 
NCASE=0 

DO 112 Te1,4 
DO 113 IP=1,5 
DO 114 Ir 1,3 
P=PVECT (IP) 
R=RVECT (IR) 
NOBS=NVECT (IN) 
NCASE=NCASE+1 
NCP=2*NCASE 
NCM=NCP~1 
ISEED1=ISD (NCP) 
ISEED2=ISD (NCM) 


CALL GNRNB (NOBS ,OBS,P,R, ISEED2 , ISDEED1) 


OBMX=-99999 .D0 
DO 321 LK=1,NOBS 


IF (OBS (LK) .GT.OBMX) OBMX=OBS (LK) 


CONTINUE 
OBLM=OBMX+10.D0 
NVAR=0 
OBMNFOBMX+5 . DO 
DO 370 MI=1,NOBS 


IF (OBS (MI) .LT.OBMN) OBMN=OBS (MZ) 


CONTINUE 
ICNT=0 
DO 380 MJ=1,NOBS 


IF (OBS (MJ) .NE.OBMN) GO TO 380 


ICNT=ICNT+1 


DAT (NVAR) =OBMN 


FRQ (NVAR) =DFLOAT (ICNT) 


IF (OBMN. BQ. OBMX) 


GO TO 395 
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395 
889 


301 
302 


543 


544 


445 


71 


72 


10 


=). 


145 


GO TO 390 
WRITE (6,889) 
FORMAT (//,10X,89('*'),/,40X,'DATA GENERATED IS GIVEN BELOW',/,39X, 
Poa) 
WRITE (6,301) (DAT(LO) ,LO=1,NVAR) 
FORMAT (20X, 'VALUES',/,20X,6('-') ,/, (5X,20 (1X,F5.0) ,/) ) 
WRITE (6,302) (FRQ(LO) ,LO=1,NVAR) 
format920X, 'FREQUENCIES',/,20X,11('-") ,/, (5X,20(1X,F5.0) ,/)) 
WRITE (6, 543) 
FORMAT (//,10X,'DATA WRITTEN ABOVE IS BASED ON FOLLOWING PARAMETERS 


1',/,10X,52 ("=") ) 


WRITE (6,544) NCASE,P,R,NOBS,NVAR 
FORMAT (/,10X, 'CASE NO.=',15.5X, 'P=',F7.4,5X, 'R=',F7.2,5X, 'N=',Ir,5 


1X, 'NO.OF VALUES=' ,13,/) 


SUM=0 . ODO 

SUM1=0. ODO 

DO 445 I=1,NVAR 

SUM=SUM+FRQ (I) *DAT (T) 
SUM1=SUML+F RQ (I) *DAT (I) **2.0D0 
XBAR=SUM/DF'LOAT (NOBS ) 

SSQR= (SUM1-DFLOAT (NOBS) *XBAR**2.0D0) /DFLOAT (NOBS-1) 
IF (XBAR.GE.SSQR) GO TO 71 
P=XBAR/SSOR 

R=XBAR**2.0D0/ (SSQR-XBAR) 

GO TO 72 

R=5.0D0 

P=0.9D0 

GO TO 114 

Q=1.0D0-P 

K=1 

TEST=Q** (K+1) / (DFLOAT (K+2) ) 

IF (TEST.LE.EPS) GO TO 8 

K=K+1 

GO TO 7 

M=K 

DETLAM+0 . ODO 

DO 10 l= 1,M 

RATION=1. ODO 

DO 12 JJ=1,N 
RATION=RATION*DFLOAT (JJ) / (DFLOAT (JJ) +R) 
DETLAM=DETLAM+O* *N*RATION/DFLOAT (N+1) 
CONTINUE 

I=1 

CALL SVS (FCN, TNOT) 

EFF 1=-FCN (TNOT) 

RSTAR=R 

PSTAR=P 

SUM2=0. ODO 

DO 145 J=1,NVAR 

SUM2=SUM2+FROQ (J) (TNOT**DAT (J) 
TXBAR=SUM2/DFLOAT (NOBS) 

CALL COMPR (TNOT, TXBAR, XBAR, RSTAR) 
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77 


88 


Die 


98 


141 


94 


29 


114 
113 
112 


P=R/ (R+XBAR) 

Q=1.0D0-P 

<=] 

TEST=Q* * (K+1) / (P*DFLOAT (K+2) ) 
IF (TEST.LE.EPS) GO TO 88 
K=K+1 


DO 51 N=1,M 

RATION=1 . ODO 

DO 52 JJ=1,N 

RATION=RATION*DFLOAT (JJ) / (DF LOAT (JJ) +R) 

DETLAM=DETLAM+)* *N*RATION/DFLOAT (N+1) 

EFF 2=-FCN (TNOT) 

IF (L.NE.1) GO TO 141 

UNOT=-DLOG (DABS (TNOT) ) 

WRITE (6,98) PSTAR,RSTAR,P,R,EFF1,EFF2,UNOT 

FORMAT (/,11X,57('&') ,/,15X, "VALUES OBTAINED AT THE END OF FIRST IT 
LERATION',/,12X, 'PTHIL', 3X, 'RTHIL' ,3X, 'PSTAR' , 3X, 'RSTAR', 3X, 'ETHIL' 
Zeros, GON 7 LX, ed, 2h, f 123, 2h,2 9.57 24,217.3,2K,89.35,5 
SNe ee 4/2, 7 & ),/) 

IF (DABS (P-PSTAR) .GT.0.001) GO TO 29 
IF (DABS (R-RSTAR) .GT.0.005) GO TO 29 
IF (DABS (EFF1-EFF2) .GT.0.01) GO TO 29 

UNOT=—-DLOG (DABS (TNOT) ) 

WRITE (6,94) XBAR,SSOR,P,R,EFF2,UNOT,L 

FORMAT (15X, 'USING THE ABOVE DATA FOLLOWING INFORMATION IS OBTAINED 
1/,/,15X,54('='),/,10X, 'SPLE.MEAN"', 2X, 'SPLE.VAR', 3X, 'EST.CF P' ,3X,' 
2ZEST.OF R',6X, 'EFF',6X, 'OP.UNOT", 3X, 'NO.OF ITR',/,10X,7(9('—") , 2X), 
3/,10X,F9.3,2X,F9.3,2X,F9.3,2X,F9.3,2X,F9.3,2X,F9.4,5X,13,//) 

GO TO 114 

I=L+1 

GO TO 95 

CONTINUE 

CONTINUE 

CONTINUE 

STOP 

END 
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SUBROUTINE GNRNB(N,OBS,P,R,ISEED2, ISEED1) 
REAL*8 P,R,OBS 

DIMENSION OBS (200) ,GAM(200) ,K (1) 
P1=SNGL (P) 

RL=SNGL (R) 

CALL GAMA (RL, ISEED1,GAM,N) 

DO 1 I=1,N 
RLAM=GAM (I) * (1.-P1) /P1 

CALL GGPOSH (RLAM, ISEED2,1,K, IER) 
OBS (I) =DBLE (FLOAT (K (1) ) ) 

RETURN 

END 


SUBROUTINE COMPR(T,TX,XB,RST) 
IMPLICIT REAL*8 (A~H,O-Z) 
COMMON P,R,DETLAM 


B=R1+XB* (1. 0D0~-T) 
TY=DABS (RL/B) **R1 
F=TY-TX 
FP1=TY 
FP2= (XB* (1. ODO-T) /B) +DLOG (DABS (RL/B) ) 
FP=FP1*FP2 
R2=R1-F/FP 
TF (DABS (R2—-R1L) .LE.0.0001D0) GO TO 66 
R1=R2 
GO TO 30 
R=R2 
RETURN 
END 
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